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We investigate the phase diagram of the half-filled SU{N) Hubbard-Heisenberg model with hop- 
ping t, exchange J and Hubbard U, on a two-dimensional square lattice. In the large-N limit, and as 
a function of decreasing values of t/J, the model shows a transition from a d-density wave state to a 
spin dimerized insulator. A similar behavior is observed at A^ = 6 whereas at A = 2 a spin density 
wave insulating ground state is stabilized. The A = 4 model, has a d-density wave ground state at 
large values of t/J which as a function of decreasing values of t/J becomes unstable to an insulating 
state with no apparent lattice and spin broken symmetries. In this state, the staggered spin-spin 
correlations decay as a power-law, resulting in gapless spin excitations at q*= (7r,7r). Furthermore, 
low lying spin modes with small spectral weight are apparent around the wave vectors q = (0, tt) 
and q = (tt, 0). This gapless spin liquid state is equally found in the SU{4:) Heisenberg {U/t ^ oo ) 
model in the self-adjoint antisymmetric representation. An interpretation of this state in terms of a 
TT-flux phase is offered. Our results stem from projective (T — 0) quantum Monte-Carlo simulations 
on lattice sizes ranging up to 24 x 24. 

PACS numbers: 71.27.-|-a, 71.10.-w, 71.10.Fd 
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I. INTRODUCTION 

SU{N) symmetric models of correlated electron sys- 
tems have attracted considerable interests in the past 
decades. For instance, those models are relevant for the 
understanding of Mott insulators with orbital degeneracy 
as described by the Kugel-Khomskii Hamiltonian Ql . For 
two-fold orbital degeneracy and at a point where the or- 
bital and spin degrees of freedom play a very symmetric 
role, this model maps onto an S'C/(4) symmetric Hub- 
bard, or Heisenbergmodel with fundamental represen- 
tation on each site Q]- It has also recently been argued 
that realizations of SU{N) Hubbard models are at reach 
in the context of optical lattices Q . 

SU{N) generalizations of SU{2) lattice fermion mod- 
els can be solved exactly in the large-N limit. Systematic 
corrections in terms of Gaussian fluctuations around the 
mean-field or saddle point solution may be computed. 
The simplifications which occur in the large- A^ limit, 
namely the suppression of quantum fluctuations have im- 
portant consequences for auxiliary flcld quantum Monte 
Carlo (QMC) simulations. As a function of growing val- 
ues of A'^ the negative sign problem inherent to stochastic 
methods is reduced thus rendering simulations more and 
more tractable. In fact, some generalizations of Hubbard 
models lead to the absence of sign problems for speciflc 
values of N and irrespective of doping |j, ^. However, 
the extrapolation from the soluble large-N limit to the 
physical N = 2 case is by no means unambiguous since 
phase transitions can occur as a function of A^. 

In this article, we will primarily concentrate on the 
half-filled Hubbard-Heisenberg model on a square lattice 
and map out it's phase diagram as a function of A^ and 
coupling strength. At this band filling, the sign problem 
is absent for even values of A^. Hence, ground state prop- 
erties can be investigated on lattice sizes ranging up to 
24 X 24 unit cells. We will show the existence of d-density 



wave (DDW) states down to A^ = 4 and of spin-dimerized 
states at A^ = 6. The most intriguing result is a possible 
realization of a gapless spin-liquid phase for the A^ = 4 
model in the Heisenberg limit. 

The SU{N) symmetric Hubbard-Heisenberg model we 
consider reads: 

H = Ht + Hu + Hj with 
Ht^-tY^ clcT + H.c. 
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,cl ) is an A^-flavored spinor, 

D:^j — clc:! and p corresponds to the band-filling. At 
N = 2, the operator identity 
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holds. Here, the fermionic representation of the spin 1/2 
operator reads S = \^s g' ^\^s,s'Cs' where a are the 
Pauli spin matrices. Thus, at A^ = 2 the model reduces 
to the standard Hubbard-Heisenberg model. 

In the strong coupling limit, U/t — > oo, and at integer 
values of pN/2, charge fluctuations are suppressed. The 
model maps onto the SU{N) Heisenberg Hamiltonian 
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the generators of SU{N) satisfying the commutation re- 
lations: 



^a,f3,i^ "^7,5,/ — ^ij \^a,S,i^7^l3 5'^,/3,i^a,i5 j 
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The representation of the SU{N) group is determined by 
the local constraint 
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N 



(6) 



In the terminology of Young tableaux the above leads to 
a tableau with pN/2 rows and a single column. In par- 
ticular, at A^ = 4, and p = 1/2 (quarter band-filling) the 
model maps onto the SU{A) symmetric Kugel-Khomskii 
Hamiltonian with fundamental representation of S'C/(4) 
on each lattice site. A study of large-N Heisenberg mod- 
els in various representations may be found in Ref. |^. 

SU{N) Heisenberg models have been considered nu- 
merically in Refs. [3,|g. Those models differ substan- 
tially from ours in the choice of the representation. On 
one sublattice the fundamental representation (Young 
tableau with one row and a single column) is consid- 
ered and on the other the adjoint representation (Young 
tableau with N — 1 rows and a single column). Based 
on Green function Monte-Carlo methods, it has been ar- 
gued this 5'[/(4) model has a spin-hquid ground state. 
However, simulations on larger lattice sizes with the loop 
algorithm have shown that the model has a broken sym- 
metry ground state |3|- In contrast, our results for the 
S'C/(4) Heisenberg model, at p = 1 and in the self-adjoint 
representation (see Eq. ((BJ) point towards an insulating 
state with no broken symmetries. 

The article is organized as follows. In the next section 
we formulate the the partition function of the model as 
a path integral over bosonic field. This formulation con- 
stitutes the starting point for both the saddle point ap- 
proximation and the auxiliary field quantum Monte Carlo 
simulations. In Section Uni we present the phase diagram 
of the half-filled model as a function of N and coupling 
constants. Finally, we summarize and draw conclusions. 



Here mAr = /3 and we have omitted the systematic er- 
ror of oder Ar^. Using the Hubbard Stratonovich (HS) 
transformation, we introduce bosonic fields to decouple 
the two body interaction terms. Let us start with the 
Heisenberg term which we write - replacing the sum over 
nearest neighbors {i,j) by a sum over bonds 6 - in terms 
of perfect squares 
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We can now apply the standard HS transformation to 
obtain: 
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In the above, we introduce a complex variable per bond: 
^b — {lb + irjb) / V^NJAt . Following the same steps, we 
decouple the Hubbard term as: 

(10) 

Using the above transformations, the partition func- 
tion in the limit At — > is given by a functional integral 
over the space and imaginary time dependent HS fields: 

Z(x fYi ^H^) n DRez6(T)DImzb(T)e-^^«*>'{">'{^». 

? b 

(11) 

The action reads: 

5({$},{z},{z-}) = /'drJ^|z,(r)p + C/^<i>H(r)/4 
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II. LARGE-N LIMIT AND QUANTUM 
MONTE-CARLO SIMULATION. 



Both the saddle point approximation as well as the 
auxiliary field QMC rely on a path integral formulation of 
the partition function. Using the Trotter decomposition, 
we write the partition function as: 



Z =. Tr [e-P"] = Tr 
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Notice that in the above definition of h{T), the creation 
and annihilation operators are not N-flavored spinors but 
correspond to spinless fermion operators. 



A. The Saddle point. 

In the large-N limit, the saddle point approximation 
SS SS 6S 
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becomes exact. Assuming time independent fields, we 
derive the mean- field equations: 
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FIG. 1: (a) The four site unit cell with lattice vectors, 
ai = 2ax and ag = 2ay and corresponding fields. Here, a^ 
and Ey correspond to the lattice vectors of the underlying 
square lattice, (b) Mean-field order parameters (see Eq. I17II 
as obtained from the saddle point equations. 

The above saddle point has been considered by Affleck 
and Marston [^ [l3|. At half-band filling, p — 1, and 
large values of t/J a d-density wave state is realized. 
This state becomes unstable towards dimerization as the 
coupling t/J is reduced. Here, we have solved the mean- 
field equations for a four site unit-cell (see Fig. QJ, thus 
allowing more freedom in the dimerization pattern than 
in llO|. With 

Z, = r,e'^■/^ (16) 

the solutions we find arc characterized by: 

n = 7-2 = rs == 7-4 ^ r^, r5 = re = 7-7 ^ r^ ^ tb 

4>1 ^ 4'2 ^ 4>3 ^ 4>4 ^ 4'A, 05 = 06 = 07 = 08 = 0B-(17) 

The values of the oder parameters are plotted in Fig. 
n As apparent the d-density wave state with r^ — tb 
and (j)A — 0B is unstable towards box dimerization below 
tc/J — 0.17 ^J. The DDW state is a semi-metal since 
gapless single particle excitations are present at wave vec- 
tors (±7r/2a, ±7r/2a). Dimerization opens a quasiparticle 



gap at all wave vectors. Hence the transition from the 
DDW to the dimerized phase corresponds to a semi-metal 
to insulator transition as shown in Fig. |21 
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FIG. 2: Lowest lying single particle excitation at half-band 
filling. The transition from the DDW state to dimerized state 
corresponds to a semi-metal to insulator transition. 

We note that the results of Affleck and Marston [13 
may be recovered by imposing: 



Z5 — Z3, Zq — Zi, Zy — Z4, and Zg — Z2. 



B. The Monte Carlo simulation. 
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The Monte-Carlo approach relies on the same formu- 
lation of the partition function. Before discussing de- 
tails of the implementation let us concentrate on our 
primary concern, namely the sign problem. In general, 
^-NS(<p,z,z) jg j^Q-j- ^ positive quantity and hence may not 
be interpreted as an unnormalized probability distribu- 
tion from which we sample field configurations. Hence, 
in the Monte Carlo method, we consider the probability 
distribution: 
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JV[(l),z,z]\e-NS{4>,z,z)\ 
and estimate the expectation value of an observable with 
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In the above, O(0, z, z) is the expectation value of 
the observable for a given configuration of fields, and 
^-NS(<t:,z,z) ^ \^^-NS(4>.z,z)\^^i5(4>.z.z)^ The denominator 

in the above equation, corresponds to the average sign: 



^e«)p. 



In the large-N limit, where the saddle point approxi- 
mation becomes exact, the average sign is temperature 
independent and equal to unity. On the other hand, it 
is known that for the SU(2) model the average sign de- 
cays as e"^^'' where V corresponds to the volume of the 
system and A is a positive constant. Hence, we can con- 
jecture that A is a decreasing function of N . This has for 



consequence that at a given temperature the sign prob- 
lem becomes less and less severe as function of growing 
values of N. We have checked this numerically for the 
quarter filled, p = 1/2, model. Unfortunately the sign 
problem for the S'C/(4) quarter-filled model - the S'C/(4) 
symmetric Kugel-Khomskii model - was still to severe 
to study the nature of the Mott insulating phase in the 
strong coupling limit. 
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FIG. 3: Average sign as a function of A'' for the quarter-filled 
SU{N) Hubbard model. The solid line corresponds to a fit to 
the form: ae~^" . 

At half-band filling, p = 1, and even values of N , 
particle-hole symmetry leads to the absence of a minus- 
sign problem. At this filling, and under the canonical 
transformation 
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the Hamiltonian h(T) transforms as: 



h{T) ^ h{T) 



such that 



Tr 



Tr 



J^^-!idTh{r) 



Hence, the above quantity is real and 
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is positive for even values of N. 

We now summarize the technicalities required to carry 
efficient simulations. Since we are interested in ground 
state properties, it is more efficient to adopt a projective 
method based on the equation: 

The trial wave function |^t) is required to be non- 
orthogonal to the ground state and /? corresponds to a 
projection parameter. For the trial wave function we 
choose the form: 



where I^t)^ is the ground state of the single parti- 
cle Hamiltonian —ty^r^x cl c^ + H.c. in the flavor a 

Hilbert space. With this choice of trial wave function, the 
action within the projection formalism takes the form: 

5 = 5o-ln(xT|re-.''o'd-''M|XT) (27) 

where \xt) is the ground state of the spinless fernuon 
Hamiltonian: — iV],^ a clc^ -I- H.c. In the simulations 
we will present in the next section, we have typically 
used PJ ~ 40 which we found to be sufficient to filter 
out the ground state from the trial wave function within 
statistical uncertainty. 

We use a finite imaginary time-step At which we have 
set to At J = 0.1. This introduces a systematic error of 
the order Ar^. Given this systematic error, it is much 
more efficient to use an approximate discrete HS trans- 
formation to decouple the perfect square term: 
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where the fields rj and 7 take the values: 

7(±1) = 1 + \/6/3, 7(±2) = 1 - V6/3 
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This transformation is not exact and produces an over- 
all systematic error proportional to (ArA)"^ in the Monte 
Carlo estimate of an observable. However, since we al- 
ready have a systematic error proportional to Ar^ from 
the Trotter decomposition, the transformation is as good 
as exact. It has the great advantage of being discrete 
thus allowing efficient sampling. 



C. The Heisenberg limit. 

We conclude this technical part with some comments 
concerning the numerical simulations of the Heisenberg 
model. At t/J — 0, Hu is a good quantum number since 
[Htuj, Hjj] — 0. Hence, in principle it suffices to choose a 
trial wave function |^t) satisfying HuI^^t) = to guar- 
antee that the imaginary time propagation converges to 
the ground state of the Heisenberg model (see Eq. |2SJ). 
On the other hand, one can relax this constraint on the 
trial wave function and implement a Gutzwiller projec- 
tion onto the Hilbert space with no double occupancy. 
We have found the second approach to be much more 
efficient. 

The algorithm we use here is very related to the one we 
have used in Ref. [IJ where a detailed technical section 
is provided. 



III. 



NUMERICAL RESULTS 
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Our results are summarized in the phase diagram 
shown in Fig. ^ Here, we consider the half-filled case 



as a function of N and t/J. For values of t/J > we 
set U = 0. The t — line corresponds to the Heisen- 
berg model where charge fluctuations are completely sup- 
pressed (see Sec. IIIC|) . In the large-N limit, the data 
stems for the mean-field calculation of the previous sec- 
tion. At A^ = 6, we essentially reproduce the saddle 
point result with a somewhat smaller value of tc/ J re- 
fiecting the instability of the DDW phase in favor of the 
spin-dimerized phase. Irrespective of the coupling t/J, 
the SU{2) model shows an insulating spin-density wave 
(SDW) state. The most interesting feature of the phase 
diagram occurs at A^ = 4. Apart from the DDW phase 
present at large values of t/J we find an insulating phase 
(solid circles in Fig. 0J) with no apparent broken symme- 
tries and no spin gap. We will argue that in this phase 
the antiferromagnetic spin correlations are critical lead- 
ing to gapless spin modes around the antiferromagnetic 
wave vector Q = (tt, tt). Furthermore, we will present re- 
sults showing that low lying spin modes with very small 
spectral weight are present around the q = (0, tt) and 
q= (0,7r) wave vectors. Before proceeding let us remind 
the reader that our simulations are carried out with the 
projective algorithm of Eq. [Ss] and hence reflect ground 
state properties. 



From the equal time correlation function So{q) = 
So{q, T = 0), we can establish the presence of long range 
order at a given wave vector. In this case So{q} scales 
as the volume of the system, V, the proportionality con- 
stant being the square of the order-parameter. From the 
imaginary time displaced correlation functions, we can 
compute spectral functions, 5*0(9,0;), by solving, 



Soiq,T) = - du;So{q,uj)e ™, 
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(31) 



with the use of the Maximum Entropy method. 

Information on gaps and the spectral weight of the low- 
est lying excitation, is obtained directly form the imag- 
inary time correlation function without having recourse 
to analytical continuation. 

Soio^r) = 5]|(*o|0(-g1IXn(?))Pe-^(^"(^"^-^») 
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In the above, |^o) corresponds to the normalized ground 
state. Ixniq)) are eigenstates of H with momentum q 
and |(*o|0(-g)|x„(g))| > 0. The gap Ao{q) corre- 
sponds to the energy difference between the first ex- 
cited state \xo{q)) and the ground state and 0{q) = 
-= ^^e"'^0(j). Finally, the residue Zo{q) corresponds 
to the spectral weight of the lowest lying excitation. 

To study the model, we have considered the following 
observables. Let us define the magnetization as 

OspiniJ) = J2 /(")clca with ^ /(a) = 0. (33) 



FIG. 4: Phase diagram of the half-filled (i.e. p = 
-2r V (cl c- )) Hubbard-Heisenberer model as a function of 

t/J. For t/J > we set [7 = 0. The t = line corresponds 
to the Heisenberg model where charge fluctuations are com- 
pletely suppressed (see Sec. lilCII . The symbols correspond 
to the parameters where we have carried out simulations and 
denote the following phases: A: Spin-dimerized phase, Q- 
DDW phase, D: Spin-density wave phase, and •: insulating 
phase with no broken lattice and spin symmetries and no gap 
to spin excitations (gapless spin-liquid (GSL) phase). 

To establish the above phase diagram, we have com- 
puted equal-time and time displaced correlation func- 
tions. Let 0{i) be an observable, with time displaced 
correlation function: 

Soi^- I r) = {Oil t)oO)) - m)){Om (29) 

and corresponding Fourier transform: 



5o(g,r)=5]e^*~5o(r,r). 



(30) 



For even values of A^ considered here, we choose /(a) = 
±1. Note that SU{N) symmetry leads to the identity 
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(34) 
where 5" are the generators of SU{N) (See Eq. ^). 

To detect spin-dimerization and DDW instabilities we 
consider respectively 



and 



Odimor(«) = Ospin(?)Ospin(« + a^:) (35) 



Oddw (i) = jx (i) - jy (i) (36) 



with current: 



c- - c- 1 (37) 
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FIG. 5: (a) Size scaling of the DDW equal time correla- 
tion functions at wave vector Q = (7r,7r). (b) Single par- 
ticle spectral function in the DDW phase on a 20 x 20 lat- 
tice. Here, we normalize the peak height to unity and along 
the (0, tt) to (tt, 0) line each spectrum satisfies the sum-rule 
J_ AujAik, lS) = ■Kn{k) — n/2. (c) Size scaling of the quasi- 

particle gap at fc = (7r/2, 7r/2) and k = (0, tt). The data shows 
the semi-metal character of the DDW phase. 



and an equivalent form for jy{i). Finally, we ob- 
tain information on single particle excitations by mea- 
suring time displaced Green functions: G{k,T) = 
— {Tcj:^{t)c-. (0)). From this quantity we can ex- 
tract quasiparticle gaps as well as the spectral function 
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t/J = 0.5,U/t = 0, iV = 4 




FIG. 6; Spin correlations in the DDW phase, a) Equal time 
spin structure factor, b) Size scaling of the spin gap at g = 
(0, tt) and q = (7r,7r). 
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A. The DDW phase 

We start our description of the phase diagram with 
the DDW phase. Fig. [SJi shows the finite size scaling of 
Sbi)w{Q)/L^ at the antiferromagnetic wave vector Q = 
(tt, tt). For t/J = 0.5 and both considered values of iV the 
data supports lim^^oo "SdowIQ)/.^^ > thus signaling a 
DDW ordered phase. One can confirm this point of view 
from the analysis on the single particle spectral function 
along the k = (0, tt) to fc = (tt, 0) line in the Brillouin 
zone (see Fig. ^\p). After size scaling (see Fig. ^) the 
data is consistent with the vanishing of the quasiparticle 
gap at fc = (7r/2,7r/2) thus confirming the semi-metal 
character of the DDW phase. 

Since the single particle spectrum has gapless excita- 
tions at the nodal points k = (±7r/2,±7r/2) we can ex- 
pect gapless spin excitations centered around the wave 
vectors 0*= (0,7r), q = (7r,0) a,ndq= (tt, tt), along with a 
power-law decay of the equal time spin-spin correlations. 
Fig. El confirms this. The spin gap vanishes at the above 
mentioned wave vectors (Fig. |H1d). The spectral weight, 
Zspin((7), of the low lying q = (0, tt) spin excitation is very 
small in comparison to q = (tt, tt). For the L = 20 lat- 
tice we have Zspin{q)/Sspin{(t) — 0.006 at g = (0,7r) and 
■Z^spin('7)/'5'spin('?) — 0.5 bX q = (tt, tt). The data on which 
this statement is based stems from Fig. 112b . The size 
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FIG. 7: Quasiparticle (a) and spin (b) gaps at A'^ = 6 as a 
function of coupling. For tj J > we set [/ = 0, and tj J = 
refers to the Heisenberg model (see Sec. IllCfl . 



scaling of the static spin structure factor at Q = (tt, tt) is 
consistent with a power-law decay of the staggered spin- 
spin correlation functions: iSspinlQ = ('t, 7r))/L^ oc L^^ . 
It is interesting to note that even though gapless spin 
excitations are present sX q — (0, tt) no sign of a cusp in 
the spin structure factor at this wave vector is apparent 
(See Fig. EJi). This signals a large exponent for the (0, tt) 
spin correlations. Note that in the saddle approximation 
both the (0,7r) and (tTjTt) spin correlations decay as r~^ 
and no cusp feature in the static spin structure factor is 
apparent 



B. N=6 

To study the stability of the DDW phase from tj J — 
0.5 to the Heisenberg point, it is convenient to consider 
the quasiparticle gap at fc = (7r/2,7r/2) (see Fig. EJi). 
For values of tj J < 1/5 the data is consistent with the 
opening of a quasiparticle gap at A; = (7r/2,7r/2). The 
opening of the quasiparticle gap is accompanied by the 
opening of a spin-gap {see[7jp). Hence, the data suggest 
that for values of t/J < 1/5 at iV = 6 we have entered a 
spin dimerized phase, with spin and charge gaps. 

We confirm this point of view by looking into the dimer 
correlation functions (see Fig. |Sli). Those correlations 
are dominant at wave vector q — (tt, 0) in agreement 
with the mean-field dimerization pattern. To establish 
long-range order we have to extrapolate S'dimer(''''7 0)/i'^ 
to the thermodynamic limit. For t/J < 1/5 gaps are 
present both in the charge and spin degrees of freedom. 



FIG. 8: (a) Equal time dimer correlation functions for the 
517(6) Heisenberg model, (b) Size scaling of the dimer corre- 
lation function at g = (0,7r). For t/J > we set f/ = 0, and 
t/J — refers to the Heisenberg model (see Sec. IllCt . 



Since the presence of gaps is equivalent to the localization 
of spin and charge degrees of freedom, we fit the finite 
size data to the form: a + br'^e~^/^. Adopting this form 
for the finite size corrections, the data is consistent with 
the onset of a spin-dimerized state for t/J < 1/5. 



C. N=4 

We now turn our attention to the SU{A) model. 
From the size analysis of the quasiparticle gap at fc = 
(7r/2, 7r/2) (see Fig. |2K) we can conclude that the DDW 
state is unstable for values of t/J < 1/4. In contrast to 
the iV = 6 case, no spin gap opens across the transition 
up to the Heisenberg limit (see Fig. |5|d). The correla- 
tion functions presented in Fig. 1101 in the Heisenberg 
limit show dominant spin-spin correlations. In accor- 
dance with the absence of spin gap no sign of long-range 
dimerization is apparent. 

The issue is now to establish the existence or non- 
existence of long range antiferromagnetic order. In 
Fig. ^2 we plot the spin-spin correlation at the largest 
distance, (L/2, L/2), on our L x L lattice as well as 
S'spinl? = ("■j'"'))/-^^ as a function of 1/L. As apparent, 
the extrapolation is consistent with the absence of long- 
range magnetic ordering. In particular, at the Heisen- 
berg point, where we have carried out simulations on lat- 
tices ranging up to 24 x 24, our extrapolated values read 
limL^oo 5'spin(r = (L/2,L/2)) = 0.002± 0.003 along with 



limr 



SspiniQ = {tt,7t))/L^ = 0.0008 ± 0.004. 
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FIG. 9: Size scaling of the quasiparticle and spin gaps at 
iV = 4. For tjj > we set U = 0, and t/J = refers to the 
Heisenberg model (see Sec. IIIC|l . 



The data is equally consistent with a power-law de- 
cay of the spin-spin correlations. Concentrating again on 
the Heisenberg point the dashed lines in Fig. Illlcorre- 
spond to the forms: Sspin{q = {t^ , t^)) / L'^ oc L~^-^^ and 



5spi„(f= (L/2,L/2))(xL- 



-1.12 



The difference between 



the two numerical values of the exponents gives an idea 
of their uncertainty. 

Hence the equal time data is consistent with an insu- 
lating phase with no apparent lattice or spin symmetry 
breaking and no gap in the magnetic excitations. It is 
now intriguing to investigate the spin-dynamics of this 
phase. Fig. 112b plots the dynamical spin-structure fac- 
tor at t/J = 0.1, f//i = on a 16 X 16 lattice. The 
data, shows several features. The gap at q = (tt, tt) is a 
finite size effect (see Fig. ^jfl). Taking this into account, 
the data is consistent with a gapless mode with linear 
dispersion around q= (tt, tt). This feature is clearly not 
surprising since the equal time correlation functions show 
critical behavior at this wave vector. As we follow this 
mode to q= (0, tt) the line shape becomes very broad and 
spectral weight seems to spill down to low energies. This 
is especially apparent on the intensity plot for which we 
have used a logarithmic scale (see Fig. ll2bV Since the 
spectral weight of the low- lying modes around q= (0, tt) 
is very small, it is desirable to confirm the above state- 
ment. To this aim we plot in Fig. 112b the imaginary time 
displaced correlation functions, Sspin{q,T), at q — (0,7r) 
and q = (tt, tt) both in the gapless spin-liquid phase at 
t/J = 0.1 and for comparison in the DDW phase at 
t/J = 0.5. Let us start with the DDW phase where 
we were able to show the presence of gapless spin modes 
around the (tTjTt), (0, tt) and (7r,0) points in thermody- 



FIG. 10: Equal time dimer and spin correlation functions for 
the S'f/(4) Heisenberg model. 
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FIG. 11: Size scaling of the equal time spin-spin correlation 
functions. For t/J > we set C/ = 0, and t/J — refers to 
the Heisenberg model (see Sec. IIIC|l . 



namic limit (see Fig. EJ)). On the L = 20 sized system 
considered in Fig. I12h one sees that both the q — (0, tt) 
and q = (tt, tt) correlators decay assymtotically with the 
same exponential form, thus signalling low energy spin 
excitations at both wave vectors. Of course there is a 
big difference in the prefactor mutiplying this exponen- 
tial decay. This replects the fact that the spectral weight 
of the low-lying (0, tt) excitation is much smaller than 



(a) S(q,uj), N = A,t/J = 0.l 

T 




FIG. 12: (a) Dynamical spin structure factor on a 16 x 16 
lattice. We have normalized the peak height of each spec- 
trum to unity so as to put forward the overall shape of the 
dispersion relation. The reader however has to bear in mind 
that the weight under each spectrum is equal to the equal 
structure factor Sspin{q) which is peaked at q = (tt, tt) and 
vanishes at g = (0, 0). (b) Intensity plot of the data of (a) on 
a logarithmic scale, (c) Imaginary time correlation functions. 
Here, one can see that without having recourse to analytical 
continuation, the bare imaginary time data supports the exis- 
tence of low lying modes at g*= (0, tt) (See text). At t/J — 0.5 
{t/J = 0.1) we consider a 20 x 20 (16 x 16) lattice. 



that of the q= (tt, tt) excitation. Let us now turn our at- 
tention to the gapless spin Hquid phase at t/J = 0.1. As 
apparent the data shows a very similar behavior at the 
exception that the spectral weight at g = (0, tt) is reduced 
in comparison to the DDW phase. We note that due to 
the extremely small scales involved in the q= (0, tt) data 
at t/J = 0.1 we are unable to obtain accurate results 
beyond tJ ~ 8. Hence, on the basis of this data, we 
believe that the spin-liquid phase indeed shows low ly- 
ing spin modes not only at (7r,7r), but also at (0,7r) and 
(tt, 0). Finally, we note that as we approach the Heisen- 
berg point, the spectral weight of the low lying mode at 
q= (0, tt) becomes smaller and smaller. This renders nu- 
merical detection of this feature at the Heisenberg point 
hard. 



IV. SUMMARY AND DISCUSSION. 

The half-filled SU{N) Hubbard Heisenberg model 
shows a variety of phases, which we have analyzed in 
detail. The saddle point physics - a DDW phase at large 
values of t/J which becomes unstable to a spin dimer- 
ized state - is valid down to iV = 6. On the other hand, 
the SU{2) model has a SDW insulating phase irrespec- 
tive of the coupling constant. The most intriguing aspect 
of the phase diagram, is the gapless spin liquid state in 
the SU{4:) model in the vicinity of the Heisenberg point. 
The SU{4) model has a DDW ground state at large val- 
ues of t/J. As appropriate for this semi-metallic state, 
we find gapless single particle excitations at wave vectors 
k = (±7r/2, ±7r/2). In the particle-hole channel those 
single particle excitations lead to gapless spin modes cen- 
tered around q — (7r,7r), q — (0, tt), and q — (tt, 0). Re- 
ducing the magnitude of the hopping matrix clement we 
find a semi-metal to insulator transition. In the insulat- 
ing phase the antiferromagnetic, Q = (tt, tt), spin correla- 
tions are critical and for the S'f/(4) Heisenberg model the 

data is consistent with the form S'spin(?^) ^ e^^''^\r\^^'^'^. 
This state shows no lattice broken symmetries and hence 
is a candidate for a gapless spin liquid state. In the 
particle-hole channel, the dynamical spin structure fac- 
tor points to gapless excitations at q — (tt, tt) but also 
to low-lying modes with small spectral weight centered 
around q= (0, tt) and q = (7r,0). 

It is tempting to argue that the gapless spin liquid 
phase is well described by a DDW mean-field state sup- 
plemented with a Gutzwiller projection. Requiring in- 
variance under time reversal symmetry pins the flux in 
each elementary square to tt. Clearly, the Gutzwiller pro- 
jection triggers the semi-metal to insulator transition but 
one could argue that in the particle-hole channel exci- 
tations remain gapless such that low lying spin excita- 
tions are present at wave vectors q— (tt, tt), q — (0, tt), 
and q = (tt, 0). In the SU{2) case, this variational wave 
function has been investigated in details. It turns out 
that due to a local SU{2) symmetry [lj| it is equivalent 
to a BCS d— wave Gutzwiller projected wave function 



10 



[l5J . At the particle-hole symmetric point the equal time 
spin-spin correlations of this wave function have been 
computed [l£j to obtain a power-law decay of the an- 

tiferromagnetic correlations: S{r) ~ e''^''|r|~^'^. Fur- 
thermore, the spin structure factor as computed form 
the variational wave function shows no cusp feature at 
q= (0,7r) and (f = (7r,0) 22]. This result compares fa- 
vorably with behavior of the spin-spin correlations in the 
SU{4:) Heisenberg model discussed in the present work. 
Alternatively we can ask the question of whether or 
not our results for the 5f7(4) Heisenberg model are un- 
derstandable from the perspective of the large-N saddle 
point. In the large- A^ limit and at the Heisenberg point 
one can stabilize the tt-Aux phase by adding biquadratic 
terms to the Hamiltonian 10] . Furthermore, Hemerle et 
al. 113 have argued in favor of the stability of the tt-Aux 
phase to gauge fluctuations arising from the constraint 
of no double occupancy. At the mean-field level, the tt- 
flux phase shows antiferromagnetic spin-spin correlations 

which decay as S{f) ~ e*'^''|r|^^. Including gauge fluc- 
tuations in a 1/A^ approximation reduces the mean field 
exponent [ig- On the other hand, the (0,7r) spin-spin 
correlations remain unaffected by gauge fluctuations [Tg ■ 
Those results compare favorably with our calculations at 
A^ = 4. 

To investigate the nature of the spinless liquid state 
we find in the SU{4:) Heisenberg model, it is desirable to 
investigate it's behavior under perturbations. Following 
the variational work of [I^ [l^i the particle- hole sym- 
metric point is unstable towards a Z2 spin liquid as re- 



alized for example in quantum dimer models |2G| . This 
instability is triggered by the inclusion of a next nearest 
neighbor hopping matrix element in the BCS Slater de- 
terminant. In the spin model this translates in the inclu- 
sion of a frustrating exchange coupling. Unfortunately, 
this is not accessible to the quantum Monte Carlo ap- 
proach since frustration leads to a minus sign problem. 
Another perturbation, which is accessible to the Monte 
Carlo approach, is the inclusion of a uniform magnetic 
field. Based on the mean-field description of the 7r-fiux 
phase, we can speculate that the point-like Fermi surface 
at zero field evolves to rings around centered around the 
zero field nodes. Since we are at a particle-hole symmet- 
ric point and that this symmetry is not broken under the 
inclusion of a magnetic field, the finite field Fermi sur- 
face is unstable towards magnetic ordering. Very much 
as in discussed in Ref. ^^ this produces a field induced 
transition to an magnetically ordered state. 
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